home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / ADDMAP.DOC next >
Text File  |  1988-10-15  |  9KB  |  280 lines
  1.         AN EXAMPLE OF HOW TO ADD A NEW MAP
  2.  
  3.     The following two functions bob() and init_bob() were added 
  4. to YMAPS.C to add the map called by the name "rr".
  5.  
  6. bob() /* this routine defines the map */
  7. {
  8.     int rand();/* a Microsoft random integer generator */
  9.     double temporary, c, s, angle;
  10.  
  11.     if ( rand() < C2*32768.)/* rand() returns a random integer 
  12.         between 0 and 32768, so C2 should be a number 
  13.         between 0 and 1, so that the first process will be 
  14.         carried out C2 of the time */
  15.     {
  16.         y[0] = 2. -y[0];
  17.         y[1] =    -y[1];
  18.     }
  19.     else
  20.     {
  21.         angle = C1*twopi/360;
  22.         c = cos(angle);
  23.         s = sin(angle);
  24.         temporary =    rho*(c*y[0] + s*y[1]);
  25.         y[1] =         rho*(-s*y[0] + c*y[1]);
  26.         y[0] = temporary;
  27.     }
  28. }
  29.  
  30. init_bob()/* this process initializes the constants and with the first 
  31.     line tells the program the name of the routine that defines the 
  32.     map; since C1, C2, and rho are used in bob(), htey must be 
  33.     initiallized in this routine; if they are not, the program will 
  34.     not let be accessed interactively and they will be given default 
  35.     values, namely -9999. */
  36. {
  37.     map = bob;/* this tells the program the name of the routine that 
  38.         is evaluated; map is a pointer to routines and bob() 
  39.         must be defined in the file prior to this line; */
  40.     rho = .65;
  41.     C1 = 135.;
  42.     C2 = .5;
  43.     /* the following 4 determine the x and y scale */
  44.     X_lower = -3.;
  45.     X_upper =  5.;     /* x scale */
  46.     Y_lower = -3.;
  47.     Y_upper =  3.;     /* y scale */
  48. }
  49.  
  50.  
  51.     The following expressions were added to YMAPMENU.C into the 
  52. routine mapmenu(). The name "rr" is name of the routine that you will
  53. use to access thses routines. Use lower case letters.
  54.     TESTMAP2 is used for expressions that will appear in two ways: 
  55. the menu and the header. The header appears in the Main Menu and the 
  56. Parameter Menu and the printout. The text (in quotes) includes the map 
  57. name so the user knows what to call the process. This first instruction 
  58. results in the following line:
  59.  
  60. RR  at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0
  61.  
  62. TESTMAP (no "2") is for material that appears only in the header. 
  63. Note that TESTMAP is used twice since two lines of header material 
  64. is to be used. THe result is the following header:
  65.  
  66. RR  at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0
  67.     The probability of a flip (x -> 2-x, y-> -y) is C2
  68.     Otherwise rotate C1 degrees about 0,0 and multiply x and y by rho (<1)
  69.  
  70. The source code for the program is:
  71.  
  72.  
  73.     TESTMAP2("rr")        fputs(
  74. "RR  at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0\n"
  75.         ,output);
  76.  
  77.     TESTMAP("rr")
  78.                  fputs(
  79. "    The probability of a flip (x -> 2-x, y-> -y) is C2\n"
  80.             ,output);
  81.     TESTMAP("rr")
  82.                  fputs(
  83. "    Otherwise rotate C1 degrees about 0,0 and multiply x and y by rho (<1)\n"
  84.             ,output);
  85.  
  86.  
  87.  
  88.  
  89. The following was added to YMAPMENU.C into the routine init_map() so that 
  90. the map could be initiallized (using the initialization routine which I 
  91. called init_bob()); initialization is important in part because it tells 
  92. the program what the name of the function is (bob in this case);
  93.  
  94.      TESTMAP("rr")
  95.     {
  96.         init_bob();
  97.         continue; /* causes exit from switch */
  98.     }
  99.  
  100.  
  101.  
  102.     THE VECTOR y[] COMPUTING LYAPUNOV EXPONENTS
  103.  
  104.     The vector y[] can have dimension up to EQUATIONS which is now 21. 
  105. The vector y includes the phase space variables and time for differential 
  106. equations and all the coordinates of the tangent vectors used in computing 
  107. Lyapunov exponents. 
  108. 1: the dimension of the phase space (not including time); this is also 
  109.         the number of coordinates of the Lyapunov vectors; 
  110. 2: the number of Lyapunov vectors to be computed initially; this can be 
  111.         changed while running;
  112. 3: lyapzero is the number of coordinates exclusive of the lyapunov vector 
  113.         coordinates; the first coordinate of the first Lyapunov 
  114.         vector is y[lyapzero];
  115.  
  116. The following are used in the routine initHenon in the file YHENON.C 
  117. to specify these numbers.
  118.  
  119.     zeroth = 0;    /* the first coordinate that is a phase space 
  120.         coordinate; it is assumed that the "vec_dim" phase space 
  121.         coordinates are together so that coordinate 
  122.             zeroth + vec_dim 
  123.         is the first coordinate following the phase space coordinates;
  124.         the non phase space coordinates may either correspond to 
  125.         the time variable or the time modulo some number such as 
  126.         twopi (which is defined) or may be some other variable of 
  127.         interset such as the total sum of the values of some 
  128.         coordinate;
  129.     vec_dim = 2;   /* the dimension of the 
  130.         Lyapunov vectors = phase space dim; this number is used in a 
  131.         number of places; it is used in computing lyapunov exponents; 
  132.         it tells what the dimension of the Lyapunov vector is; 
  133.         it is similarly used in Newton method calculations; it is also 
  134.         used when the distance between two vectors is needed; the 
  135.         distance usually useds only the phase space coordinates */
  136.     zeroth = 0;    /* the first coordinate that is a phase space 
  137.         coordinate; it is assumed that the "vec_dim" phase space 
  138.         coordinates are together so that coordinate 
  139.             zeroth + vec_dim 
  140.         is the first coordinate following the phase space coordinates;
  141.         the non phase space coordinates may either correspond to 
  142.         the time variable or the time modulo some number such as 
  143.         twopi (which is defined) or may be some other variable of 
  144.         interset such as the total sum of the values of some 
  145.         coordinate;
  146.     num_lyap = 0;  /* This is set only if you give the program the 
  147.         equations for computing Lyapunov exponents;
  148.         the number of Lyapunov numbers 
  149.         to be computed <=vec_dim; if this is not set; you will not be 
  150.         able to change its value later from within the program; if 
  151.         you do not tell the program what the partial derivatives are 
  152.         then you do not set num_lyap; the fact that it is not set to 
  153.         any value means it takes the default value of -9999 */
  154.     lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
  155.         zeroth lyapunov vector */
  156.  
  157. Computation of Lyapunov exponents is done in such a way that very 
  158. little is wasted when Lyapunov exponets are not wanted. The following 
  159. is a variant of the actual routine Henon() in file YHENON.C.
  160.  
  161.  
  162. Henon()
  163. {
  164.     int i;
  165.     int base;
  166.     double y_temp[6];
  167.  
  168.     y_temp[0] = rho -y[0]*y[0] + C1*y[1];
  169.     y_temp[1] = y[0];
  170.  
  171.     /* the following "for(...){...}" is for computing tangent vectors; */
  172.  
  173.     for (i = 0; i < num_lyap; i++)/* if the number of Lyapunov 
  174.         exponents to be computed num_lyap  is currently 0, 
  175.         then this part of the routine is skipped */
  176.     {
  177.       base =  lyapzero + vec_dim*i;
  178.       y_temp[base] 
  179.         =  -2*y[0]*y[ base ]
  180.              + C1*y[ base + 1];
  181.       y_temp[base+1] 
  182.         =  y[ base ];
  183.     }
  184.     for (i = 0; i < lyapzero + vec_dim*num_lyap; i++)
  185.       y[i] = y_temp[i] ;
  186. }
  187.  
  188.  
  189.  
  190.                   ADDING A DIFFERENTIAL EQUATION
  191.  
  192.     The program knows that a differential equation is going to be 
  193. used when it finds the differential equation step size has been set, 
  194. that is the variable "step",e.g.:    step = 1./100.;
  195.  
  196.     For differential equations, the program must be told in the 
  197. initialization file what the name of the vector field is, DELorenz in 
  198. this case. Thus the function pointer DE is set in the Lorenz 
  199. initialization with the line:
  200.     DE = DELorenz;  
  201. The pointer DE is called by rungekutta(), or any other solver used. 
  202. The Lorenz vector field is defined as follows; Y[] denotes the 
  203. position at which the vector field will be evaluated and k[] is the 
  204. vector field. Notice that there are 4 coordinates, one of which is time.
  205. Just as with y[], k[] must have coordinates for phase space, time, and 
  206. all lyapunov vectors.
  207.  
  208. DELorenz(k,Y) 
  209.     double Y[],k[];
  210. {
  211.     int i,base;
  212.     double x1,y1,z1;
  213.     double J00,J01,J10,J11,J12,J20,J21,J22;
  214.  
  215.     x1 = Y[0];/* this and the following 2 lines are for speeding
  216.         computation ever so slightly and for increasing 
  217.         readability */
  218.     y1 = Y[1];
  219.     z1 = Y[2];
  220.  
  221.     k[0] = -sigma*(x1 -y1);
  222.     k[1] = rho*x1 -y1 -x1*z1;
  223.     k[2] = x1*y1-beta*z1;
  224.     k[3] = 1; /* time */
  225.  
  226.     if (num_lyap > 0)  /* the following are the partial derivatives 
  227.         of the vector field */
  228.     {
  229.         J00 = -sigma;
  230.         J01 =  sigma;
  231.         J10 = rho-z1;
  232.         J11 = -1.;
  233.         J12 = -x1;
  234.         J20 = y1;
  235.         J21 = x1;
  236.         J22 = -beta;
  237.      }
  238.     for (i = 0; i < num_lyap; i++)
  239.     {
  240.       base =  lyapzero + vec_dim*i;
  241.       k[base + 0] =
  242.         J00*Y[base] + J01*Y[base+1];
  243.       k[base + 1] =
  244.         J10*Y[base] + J11*Y[base+1] + J12*Y[base+2];
  245.       k[base + 2] =
  246.         J20*Y[base] + J21*Y[base+1] + J22*Y[base+2];
  247.     }
  248. }
  249.  
  250.  
  251. initLorenz()
  252. {
  253.     int DELorenz();/* this is not needed if DELorenz() is defined 
  254.         in the same file as this routine and is earlier in the 
  255.         file than this routine. */
  256.             
  257.  
  258.     vec_dim = 3;   /*the dimension of the 
  259.         Lyapunov vectors = phase space dim*/
  260.     num_lyap = 0; /*the number of Lyapunov numbers 
  261.         to be computed <= vec_dim  */
  262.     lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
  263.         zeroth lyapunov vector */
  264.  
  265.     X_upper =  25.0;     /* x scale */
  266.     X_lower = -25.;
  267.     Y_upper =  60.;     /* y scale */
  268.     Y_lower =  0.;
  269.  
  270.     Y_coord = 2;
  271.  
  272.     sigma = 10;
  273.     rho = 28.0;
  274.     beta = 8./3.;
  275.  
  276.     step = 1./100.;
  277.     DE = DELorenz;   /* DE is a pointer to a function   */
  278. }
  279.  
  280.